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Abstract 

We study historical calibration of one- and two-factor models that are known to describe relatively well 
the dynamics of energy underlyings such as spot and index natural gas or oil prices at different physical 
locations or regional power prices. We take into account uneven frequency of data due to weekends, holidays, 
and possible missing data. We study the case when several one- and two-factor models are used in the joint 
model with correlated model factors and present examples of joint calibration for daily natural gas prices 
at several locations in the US and for regional hourly power prices. 
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1 Overview 

As is well-known the calibration of valuation models is a crucial step towards the realistic valuation of financial 
contracts. In this paper we focus on historical calibration of one- and two-factor models that describe the 
dynamics of energy underlyings such as spot prices of natural gas at a physical location or regional market 
power prices. These models are known to perform relatively well for energy underlyings, for example, see 
[1 or PJ. From these simple models the joint model of simultaneous evolution of several underlyings can be 
constructed by specifying the term structure of correlations between model factors. This joint model can be 
used for valuation of energy derivative contracts that depend on several underlyings, such as swing, transport, 
or structured contracts that are abundant in the energy OTC markets. These models are extensively used and 
tested by various FEA clients for valuation and risk management of OTC and standard energy contracts. 

We choose a simple robust approach for the joint calibration of these models. First, we calibrate each model 
separately based on the historical time series of prices. After that, using the calibrated parameters and the 
time series of prices, we find the time series of stochastic factors of these models and compute the correlation 
between these factors to obtain the correlation term structure of the joint model. 

As the result of calibration, we obtain volatility and correlation term structures and mean-reversion rates 
together with confidence intervals for these estimates. We also obtain historical mean-reversion level term 
structures, which include the market price of risk and are less relevant for valuation since it is done in the 
risk-neutral framework. After calibration is performed, we run the statistics analysis of model factors time 
series to see if they satisfy the initial modelling assumptions and if the model is a good fit for the data. 

One of our goals is providing rigorous calibration procedures that take into account the structure of financial 
energy data, in particular, the fact that data is not available on weekends and holidays, making the data 
frequency or, equivalently, the data step uneven. The possibility of missing data only exacerbates this effect. 
In this context, the estimation of mean- reversion rates in the models considered becomes non-trivial and we fully 
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address this issue. We also take into account seasonal effects that are typical for energy markets by allowing 
term structures of model parameters such as volatility and correlations. We provide confidence intervals for 
the estimates to see if they are statistically significant, which should be defined on a case-by-case basis. For 
some data sets, in order to increase the robustness of estimators, the granularity of the model term structures 
should be chosen in such a way that it leads to smaller confidence intervals, for example, the correlation term 
structure can be chosen to be flat and the volatility term structure can be chosen to be seasonal rather than 
monthly. These are the modelling choices that affect the value of calibrated parameters and the robustness of 
estimators. 

The structure of the paper is as follows. First, we describe the one- and two-factor models we use for 
valuation and their calibration. Then we describe how these models are joint together for the multi-asset case 
valuation and provide calibration results of natural gas prices at several location for the US market. We also 
provide an example of calibration with regional ERCOT power prices. 

Acknowledgements We would like to thank Angelo Barbieri, Tsvetan Stoyanov, and Maksim Oks at 
MSCI Research Valuation and Energy group for numerous fruitful discussions and comments. 



2 One-Factor Mean-Reversion Model 

The one-factor log-normal mean-reversion model is defined by the following stochastic differential equation 
(SDE) 

dlogS t = (8 t -alogS t )dt + a t dW u (1) 

where 9t is derived from the no-arbitrage condition Eo [St] = F(0,t) under the risk-neutral measure. 
Analogously, the one-factor normal mean-reversion model is defined by the following SDE 

dS t - (6 t - aS t )dt + a t dW t (2) 

with the same no-arbitrage condition. 

Further on, we will focus on the log-normal mean-reversion model, the results for the normal model are 
very similar and can be easily obtained in the same way. 



2.1 Risk-Neutral Valuation 

We denote X t — log S t ■ 
It is easy to get 

X t = e- a ^X s + f e- a ^9 u du+ f e- a ^a u dW u . (3) 

J s J s 

From the no-arbitrage condition we have 

f e -«(*-«)6» u rf u = logf(0,t)-e- at logF(0,0)-i / e - 2a(t -< ] a 2 u du. 
Jo 2 J Q 

It follows that 

S t = F(0, t) exp (-\v 2 {0, t) + J e- a ^a u dW u 



with 



V 2 ' 



F(t,T)=nO,T)expQ(l- e -^)) e -^)V 2 (0,t)) (^A ' 



(4) 
(5) 



2 (0,t)= / e- 2a[t - u) G 2 u du. (6) 
Jo 

If we define the evolution of the forward curve to be F(t, T) = E 4 [St], we get 



(7) 



,*) 

Thus, we see that the dynamics of the whole forward curve is explicitly defined via the dynamics of the spot. 
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2.2 Calibration 



We use historical data for calibration of the model parameters. Now we are considering the model in the real- 
world measure, hence, the mean-reversion level function 9 t contains the market price of risk and is different 
from the result obtained in the previous section [2. II 

Different granularity can be assumed for the term structures of mean-reversion level and volatility. Without 
lack of generality, let us assume that the granularity of both term structures is monthly, although, our results 
generalize to an arbitrary time indexation scheme. 

In addition, we assume yearly periodicity of volatilities, i.e. we assume that volatilities in the same calendar 
months in different year are the same. On one hand, this assumption allows us to capture seasonality effects, 
on the other hand, it enables us to use several years of data for volatility estimates making them more robust. 

Thus, the term structure of volatilities is parameterized by 12 parameters a m , with m = 1, .., 12, and the 
term structure of mean-reversion level 0( m ,y) is parametrized by a pair of parameters (m, y), with m = 1, .., 12 
and y ranging over the years specified in the data set, e.g. y — 2008, 2009, and 2010. 

Note that we are considering data sets with varying time steps since data is not available on weekends and 
holidays. 

Let us denote X t = log St for log-normal and X t — St for normal model. Solving the above SDEs explicitly, 
we obtain 

1 _ „ — adt / 1 — p — 2adt 

X t+d t = e- adt X t + 9 t +<W e t , (8) 

a V 2a 

where et is a standard normal variable. 
Let us denote 



Vt = e~ adt , 



g—adt 



«t - , (9) 



It 



1 - e 



-2adt 



2a 

The conditional probability of going from X t to X t +dt is given by 

(V W Q\ 1 f ( X t+dt - VtX t - OtKt) 2 , 

p(X t+ dt\Xt,a,a t ,9 t ) = -= exp (10) 

V27rcr t 7t \ 1<T tit 

The Maximum Likelihood Function is given by 

N-l 



C(a,a,6) = lo &P( x tk+i\ x t k ,a,(7t k ,0t k ), (11) 



k=l 

where N is the number of data in the provided time series for the underlying. 
Differentiating w.r.t. 6t m ,y)i we obtain 

„ _ St fc _ ie ( m ,y) ( X tk ~ r]tk-l X t k -i) 

V{m,y) v"-* 



(12) 



Thus, we see that the estimator of 9t depends only on the mean-reversion rate parameter. 
Differentiating w.r.t. a mi we obtain 

1 f Xt " ~ %fc-i^tfc-i ~ 6{m,y){a)K tk _ 1 \ 2 
ff ™ = ^) L { ^ ) ■ ( 13 ) 

V ' t fc -iG(m) V lt "- 1 7 

Therefore, we see that the estimator of at also depends only on the mean-reversion rate parameter. 
Thus, we have proven the following result. 

Theorem 2.1 The Maximum Likelihood estimators of mean-reversion level and volatility term structures for 
the one-factor mean-reversion model are given by 
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, , ^t h -ie(m,») \ X tk flt k -i X t k -i) , , 

l( m ,y)(a) = ^ (14) 



'-<■> - ivr E ( ^ J ' (lo) 

tfc-i£i.™j 

TTie mean-reversion rate can be found by finding the global maximum of the function 

JV-l 

= H lo gK^t*+il^t*> a > CT u(a)>M a ))- (!6) 
fc=i 

or 

C ( a ) = 2^ - lQ g(^( a )7tJ 2 2 7372 . ( 17 ) 

where parameters 7 and k are defined by {P|) and aZso depend on a. 

The maximum of this function can be found by using the Brent one-dimensional search method with the 
initial guess provided by the estimate of the mean-reversion rate based on the assumption of constant time 
steps (or, equivalently, on the assumption that the regression coefficients are slowly varying with time). 

We note that, as usual, the unbiased variance estimator for volatility should be slighlty adjusted 

/s 1 ( x t k - i]t k - 1 X tk _ 1 - 0(m,y)(a)«t fc _i \ 2 
o-m(a) = — > . (18) 

N{m) ~ 1 t k _Mm) V J 

2.2.1 Constant Time Step 

Let us see how the results can be simplified if time steps are constant. 
In this case, we have 

9(m ' v) = N(m,y) ^ Xtk ~ Mm, y) ^ 

We introduce the mean level function f t such that 

t = d t \og f t + alog f t . (19) 

Then, it follows that 

1/N(m,y) 

ft=\ n x 

\tfe-ie(m,a) 

and 

1/N{m,y) 

ft+dt = I \\ S tk 

\tk-i£(m,y) 

for t G (m,y), i.e. the mean level is the geometric average of prices in the month (m, y). 
We can rewrite the model process (p} as 

dlog{S t /ft) = -alog(S t /ft)dt + atdWt. (20) 

Denoting x t = \og(St/ ft) and integrating we obtain 

xt+dt = KX t + (Ttjet, (21) 
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where k and 7 are constant. We bucket the time series into calendar months over which at is also constant 
a m . Let us also denote N{m) by n. Then, we perform the regression to obtain 

k = nS *v- S * S v ^ (22) 



where 



Sx — ^ X tk-1> 

tfc-ie(m) 

Sy = X] Xtkl 

tfc-ie(m) 

S X x = y ' X tk-i> (^) 

tfc-lG(m) 

^I/ = Xt k-l X tk' 

tfc_i€(m) 

A = nS xx — SL 



with the standard error given by 



^ - S ra /A. (24) 
It follows that the mean-reversion rate is given by 

a = -log(/c)/dt. (25) 

When the time steps are not approximately constant, the following estimate gives a result close to the MLE 
result. The day count C of the time scries can be computed by 

C = 3657V{num of points in time series}/7V{time length of the time series}. (26) 

The estimate of the mean-reversion rate a can be obtained by following the above regression procedure. After 
that, the estimate of the mean-reversion rate should be adjusted by the day count, i.e. 

a = aC. (27) 

However, we note that for data with varying time steps this simple estimation procedure can fail for some 
time series, e.g. we encountered the cases when it fails for the time series of weekend power prices. 

2.3 Zero Mean-Reversion Rate Limit 

In the limit when the mean-reversion rate goes to zero, all the formulae derived for the one-factor mean-reversion 
model are well-defined and the one-factor mean-reversion model reduces to the Black-Scholcs model 

d log S t = ntdt + <T t dWt, (28) 

where n t = lim a ^ @t- 

The following proposition naturally follows 

Proposition 2.2 The MLE estimators of the mean-reversion level and volatility term structures in the Black- 
Scholes model are given by 

T.t k ^(m,y) ( X t k 

M(m,w) = -77— (29) 



and 



N(m) ^ dt k 
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2.4 Calibration Procedure Steps 



Here we summarize the calibration procedure. 

1. Find the estimate of mean- reversion rate do by performing regression under the assumption of the 
constant data step. (Or use an arbitrary reasonable initial estimate.) 

2. Find the maximum of the ML function using the one-dimensional search procedure with the found initial 
value of the mean-reversion ao- 

3. Compute the estimators for the term structures of the mean-reversion level and volatility. 

4. Find the error estimates for these estimators. 

5. Find the residual /mo del factor time series of the model. 

6. Compute the statistics of the residual time series to see if the model is a good fit to the historical data. 
Accept or reject the modelling hypothesis based on the Jarque-Bera and/or Kolmogorov-Smirnov tests. 

We describe the computation of the error estimates and statistics of the residuals in the following sections. 



2.5 Goodness of Model Fit 

After the calibration is done and all the parameters are found, we can find the residual time series of the model 

Xt+dt — VtXt — t K t , . 

e< = • (31) 

a tit 

and compute the first four moments of the residual time series (see section [5] for examples). We found that a 
typical situation for the real energy price data is that the first three moments are very close to normal (0,1,0), 
however, the fourth moment is different from 0, i.e. the distribution of the residuals is leptokurtic. Thus, as 
expected, the distribution of residuals for real energy data has fat tails. 

We also compute the Jarque-Bera and Kolmogorov-Smirnov statistics of the residual time series to have a 
more rigorous tests of normality of residuals and the goodness of fit of the model. 



2.5.1 Jarque-Bera statistics 

The Jarque-Bera statistics of a time series Xi with n data points is given by 



JB = ^ + \k*) (32) 



where S is the sample skewness and K is the sample kurtosis 

Q _ 

15 — ,3 ' 



^2 

K=^-3. (33) 



Here, \ik is the kth moment of the time series 

1 



Hk = ~ Y\{xi -x) k , (34) 
n z — ' 

i=l 



where x is the sample mean 

x = — y Xi. (35) 



1 ^ 



77 



= 1 



It is well-known that the JB statistics have asymptotically chi-square distribution with two degress of 
freedom, however, this approximation is good only for large sample sizes (> 2000). Therefore, for typical sizes 
of financial time series, the results of Monte Carlo simulations should be used. We used the implementation of 
the Jarque-Bera test provided in Alglib project, which is freely available online under a general GNU license. 
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2.5.2 Kolmogorov-Smirnov statistics 

We use a simple Kolmogorov-Smirnov test. First we compute the Kolmogorov-Smirnov statistics of the time 
series x\ with n data points 

D = max (\i/n-N(xi)\,\{i + l)/n-N{xi)\). (36) 

i=Q,..,n — 1 

Then, the p- value is given by, for example, see [10] 

p = Prob(L> > D ) = Q KS ((Vn + 0.12 + 0.1 l/Vn)D ) (37) 

The following Kolmogorov-Smirnov test with a confidence level naturally follows. If Dq > D a , where D a — 
— ®), then the null hypothesis that the distribution of residuals is normal, can be rejected with 1 — a 
confidence. 

3 Confidence Intervals for the Volatility, Correlation, and Mean- 
Reversion Rate Estimates 

We use well-known results for estimation of the confidence intervals for volatility and correlation estimates, for 
example, see PQ an d [5]. 

3.1 Confidence Intervals for the Volatility Estimates 

Let us simplify the notation and denote N(m) by n and a m by a. Then 

P [(Tib <(T< <T ub \ = 1 - 

07& = o\/-2^ — - — , (38) 

V X(n-l;a/2) 



Tl- 1 



(Tub = (T a 



2 

^(n-l;l-a/2) 



where a is the confidence level and X? n -i- a /2) ^ s ^ ne va l ue °f the chi-square distribution with n — 1 degrees of 
freedom and a confidence level 1 — a. 

3.2 Confidence Intervals for the Correlation Estimates 

Let us simplify the notation and denote N(m) by n and p m by p. Then, using Fisher's transformation, we 
have the following bounds on the correlation estimates 



P [Plb < P < Pub] = 1 - 



O, 



cxp (2z lb ) - 1 

Plb = Z 77; — TTT' ( 39 ) 



Pub 



exp (2zi b ) + 1 ' 
1 - exp (-2zi b ) 
1 + exp (-2zib) ' 



where 



with 



z lb = z - c/Vn-3, 

z u b = z + cj \Jn — 3 (40) 



1. fl + p 

z = — log 



2 °\l- P/ 

c = N~ 1 (l-a/2). (41) 



N is the inverse of the cumulative normal distribution and a is the confidence level. 
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3.3 Confidence Interval for the Mean-Reversion Estimate 

A standard error estimate a ao for the initial value of the mean-reversion rate ao follows from the regression 
procedure. A crude standard error estimate of the MLE estimator of the mean-reversion rate can be obtained 
in a simple way by scaling the standard error estimate coming from regression by the ratio of the MLE and 
initial value of the mean-reversion rate: 

a 

v a = — cr ao . 

a 

In order to estimate a a more rigorously, Fisher information should be employed. 

4 Two-Factor Spot-Prompt Model 

The two-factor Spot-Prompt model is a natural extension of the one-factor log-normal mean-reversion model 
which uses the prompt- month forward price (the "index") as the stochastic mean-reversion level, see [12]. The 
model is defined by the following system of SDEs 

d log S t = {9 t + a log I t -a log S t )dt + afdW t , 

dl t = I t aidB u (42) 
Covar [dB t ,dW t ] = ptdt, 

where 6 t is derived in the risk-neutral measure from the series of no-arbitrage conditions described in the next 
section. 

4.1 Risk-Neutral Valuation 

The Spot-Prompt model is a two-factor model, where the evolution of the forward curve is described by the 
rolling prompt contract. Inside each month, the expected value of the spot is given by the level of the just 
expired forward contract (the index). The spot mean- reverts to the stochastic level of the prompt contract. 
The prompt contract follows the driftless geometric Brownian motion. 

Let us denote by Tj the end of each month from the value date t = to the expiry of the contract T, i.e. 
To < < T\ < .. < T n < T < T n+ i (following NYMEX conventions, T falls on a business day three business 
days prior to the first delivery date). Let us assume that the forward curve is given by specifying the prices 
of forward contracts expiring at the end of each month F(0,Ti). We assume that the dynamics of forward 
contracts is given by the Black model: 

dF{t,T) = F(t, Ti)a F (t, Ti)dBi (43) 

for t <Ti. We make the following two assumptions that will allow us to describe the evolution of the forward 
curve using one- factor model of the rolling prompt: 

1. All forward contracts are perfectly correlated, 

2. The local volatility function for each forward contract is given by 

a F (t,Ti) = exp(-6(T, - t))a F (Ti), 

where b is a forward mean-reversion rate. Hence, we include the time-to-maturity/Samuelson effect through 
the introduction of the forward mean-reversion rate. 

Let It — F(t,Ti) for Tj_i < t < Ti denote the prompt contract at time t and <Jj(t) = a F (t,Ti) for 
Tj_i < t < Ti denote the term structure of prompt contract volatility, then, under the above assumptions, the 
dynamics of the forward curve is effectively given by the dynamics of the rolling prompt: 

dl t = Im{t)dB t . (44) 

Note that the process It is discontinuous at Ti since at every Ti the value is switched from F(Ti,Ti) at T7~ 
to F(Ti,T i+1 ) at T+. 

In the Spot-Prompt model, we assume that the following mean-reverting process describes the risk-neutral 
dynamics of the spot: 

d log S t = (0 t + a log I t -a log S t )dt + a s t dWt , (45) 
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where C o var [dB t , dWt] — Ptdt, <7s(t) is the deterministic term structure of spot volatility, and 9t is the mean 
reversion level to be derived from the following no-arbitrage assumption: 

E T JS t ] =I T - =F{T i ,T i ) (46) 

for Ti < t < Ti+i. This assumption simply means that there is no riskless profit to be made by buying index 
and selling spot or vice versa, also see pQ. 

We note that the initial no-arbitrage condition Eo [St] = F(Q,t) holds by the tower law property and the 
fact the index process is a martingale. 

In order to find the level 0(t), we derive the following equality (we omit the details of this derivation): 



e -<^6{u)du = log (M^j - e -<*-> log [S * 



Z Z Z I „ 



for Ti < s < t < Ti+i, where 



e~ a{t - u) ai{u)as{u)p{u)du (47) 



n.s./i \ •• (^); us i 



and Var s is the conditional variance at time t given the information at time s with s < t. Computing V(s,t) 
explicitly, we have 



t 

2/ 



V(s,t) = j a z {u)du (49) 
with 

a 2 {u) = <t 2 s {u) + oj(u) - 2as(u)ai(u)p(u) (50) 



We note that (j47|) follows from the dynamics of the underlyings (|44|) and (j45j) and is true for any no-arbitrage 
assumptions. From (|4T|l it follows that 

for Ti <t < Tj-|-i. From the dynamics of the rolling prompt, we compute 

E Ti [7t] = 7 Ti (52) 
for Ti<t<T i+x . Therefore, from @BJ and it follows that 

0(i) = I F ( Ti ,t)-^ (53) 

for ^ <t <T i+1 . 

From ()53|) we see that the Spot-Prompt model with the imposed no-arbitrage assumption (|46|1 can be easily 
implemented in the Monte Carlo simulation framework, see [3]. 

St is a lognormal variable. For the simulation, we use the expected value of log St 

E s [log S t ] - fU a(t - s) logS s + (1 - e-°( t - s ))log7 s - (54) 
±V(T it t) + \e-^-^V{T llS ) + ±V(s,t) - ivar s [log St] (55) 

and the variance of log S t 

Var s [log S t ] = J a 2 {u)du~ 

2/ e- Q ( t - s) (cr2( u )_ cr7 ( M ) CTs ( M )p( u ))d W + V r (s,t) (56) 
J s 
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for T{ < s < t < Tj_|_i. We note that the spot process is continuous inside each month but is discontinuous at 
Ti in general (compare with [8]). 

We also note that by the tower law and the fact that the index is a martingale it follows that 

E s [logS t ] = F(s,T i ) (57) 

for s <Ti<t < T i+ i. 
4.2 Calibration 

In the real-world measure the model can be written as 

d log S t = (6 t + a log I t -a log S t )dt + afdW t , 
d log I t = p t dt + o{ dB t , 
Covar [dB t ,dWt] = ptdt, 

We can rewrite it as 

dlog(S t /I t ) = {6t-alog(S t /It))dt + atdW u 

dloglt = mdt + a'dBt, (58) 
Covar [dB t ,dW t ] = ptdt, 

where dWt = [crfdWt — erldB t )/er t and, therefore, 

CT t = ( cr t)+ ( CT t ) - 2(J t VtPt, 

Pt = {p%Pt-oi)/<T t) (59) 
t - Ot-pt- 

It follows that effectively we have one-factor mean- reversion model on the quotient of St and It ■ Thus we 
can apply the results of the previous chapter to get the mean-reversion rate a, the term structure of 9t and a t . 

We also apply the results of the section on the Black-Scholes model and ([50]) to get the term structure 
of pt and g\. 

Thus, we have the following result. 

Proposition 4.1 The mean-reversion rate estimate for the Spot-Prompt model follows from the MLE estimate 
for the one-factor mean-reversion model that describe the dynamics of the quotient process of spot St and index 
It- 

Next, we proceed to finding the term structure correlation between the model spot and prompt factors and 
the term structure of spot volatilities. 

Denoting X t — \og(S t /It) and Y t — log(/ t ), we have 

Xt+dt = rjtXt + §K t + <T t jt£t, 

Yt +dt = Yt + Htdt + o-lVdtt,t, (60) 
Let us introduce the following time series: 

X t = {Xt+ d t-e~ adt X t -0tm)hu 

Yt = Yt +d t-Yt- ptdt. (61) 
Then, the covariance between these time series is given by 

Covar[X t , Y t ] = a\ [of p t - a*). (62) 
From here, since we know a\, we find erf pt- Then, plugging it in (|59|) . we find erf , 



af = ^at-{aly + 2afalpt. 
Finally, dividing the product erf pt by as, we get pt- 
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4.2.1 Spot Factor Time Series 



Rewriting (|60[) . we have 

X t +dt = VtXt + 0n t + o-f 7 f ef - of 7^, 

Y t+dt = Yt+fitdt + alVdi^t, (63) 

where ef , v t , and £t are correlated standard normal variables. We need to find time times series ef for the 
computation of model factor correlations when the two-factor model is included in the joint multi-asset model. 
First, we compute the correlation between v t and £t 



We notice that the function on the right is very close to identity even when the mean-reversion rates are 
very high, e.g. when adt is 0.5, the correlation is 0.99, when adt is 1, the correlation is 0.96. Thus, we can 
approximate v t by £t, to have 

4 ~ ( X t+dt ~ VtX t - 6 t K t + -¥=(Y t+dt -Y t - mdt)) . (65) 



of 7t V Vdt 

We ran tests to see that this approximation works quite well, which follows from the statistics of the residual 
time series for different data sets. 



5 The Joint Model of Several Underlyings 

Here we describe the joint model that we use for the description of the joint dynamics of several underlyings. 
Let us assume that we need to calibrate n underlyings. Without loss of generality, we assume that the first k 
underlyings are described by the Spot-Prompt model and the rest - by the log-normal mean-reversion model. 
Then, we have 

dlogSi(t) = (fl i (t) + o i logJ i (t)*i-o i logSi(t))tft + trs < (t)iWi(t), (66) 
with 5i = 1 for i = 1, .., k and 6i = for i = k + 1, .., n, and for the index 

dlogii(f) = tii{t)dt + a h {t)dBi{t), (67) 

for i = 1, .., k. There is a correlation between factors 

Govai[dWi,dWj] = p lJ {t)dt 1 

Covax[dWi, dBj] = q tl {t)dt. (68) 

We note that the chosen granularity of volatility as { , cri i and correlation pij (t) , (t) term structures as well 
as the granularity of the mean-reversion level 9i(t) and drift fii(t) term structures is a modelling assumption 
that directly affects calibration results as well as robustness of estimates. It should be chosen on a case by case 
basis depending on the data. 

In order to calibrate the joint model, we first calibrate each model separately to get a,i, the term structures 
(jSi(t)-, and ai i (t), and 6i(t) and fiiit). After that we obtain the normalized time series of model factors 
corresponding to dWi(t) and dBi{t) driving factors. Then we compute the correlations between these factors 
to get the term structures of correlations (t) and qij (t) . 

We run the Jarque-Bcra and Kolmogorov-Smirnov tests on these time series to define how good of a fit the 
model is for the provided data. In order to improve the fit, we can introduce a simple procedure for removing 
outliers. We take the distribution of the residuals of each model and throw away some percentage of outliers, 
e.g. from 1% to 5% percent in the tails. We can think of these points as points corresponding to jumps that 
were not accounted by us in the model. We saw in our examples, that this procedure improves the fit of the 
model to the data. 

We estimate the confidence intervals for the estimates in order to see if the modelling assumption on the 
granularities of term structures is good for provided data. The correlation estimates notoriously have big 
estimation errors, so the confidence intervals for correlations should be carefully checked. If the confidence 
intervals are too big and cannot be accepted the coarser granularity should be assumed on the term structure 
and the model should be recalibrated. 
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6 Calibration Examples 



We use daily closing prices of spot natural gas at several location. In the first example we use around 10 years 
of data from 01/01/1998 to 11/25/2009 and two US location STX, in the Gulf, and M3, in the Northeast, with 
the log-normal mean-reversion model. 

In the second example, we use about 2 years of data from 2/13/2008 to 11/25/2009 at two US locations 
STX and WLA together with the index data at those locations with the spot-prompt model. 

In the third example, we use regional ERCOT hourly power prices and system load from 7/1/2004 to 
2/28/2010 with the log-normal mean-reversion model. 

6.1 Log-Normal Mean- Reversion Model Daily Natural Gas Data 

Here we use around 10 years of data from 01/01/1998 to 11/25/2009 for spot natural gas prices at STX and 
M3 with the log-normal mean-reversion model. The data step is daily. We assume that the granularity of the 
term structures of mean-reversion level, volatility and correlation is monthly. 

The results of the calibration are summarized below. The mean-reversion rate estimates and their standard 
errors are given by 





STX 


M3 


mr 


38.73 


47.48 


mrc 


2.59 


2.5 



We assume the step following interpolation of the term structures. The local volatility term structures are 
given by 



The lower bounds of the confidence i 





STX 


M3 


Nov 


1.32 


1.36 


Dec 


1.10 


1.63 


Jan 


0.73 


3.05 


Feb 


1.24 


2.15 


Ma 


0.60 


1.47 


Apr 


0.46 


0.51 


May 


0.49 


0.49 


Jun 


0.56 


0.59 


Jul 


0.54 


0.63 


Aug 


0.66 


0.70 


Sep 


0.86 


0.83 


Oct 


1.07 


1.06 


vals are 




STX 


M3 


Nov 


1.23 


1.27 


Dec 


1.02 


1.52 


Jan 


0.68 


2.85 


Feb 


1.15 


2.00 


Mar 


0.56 


1.37 


Apr 


0.43 


0.47 


May 


0.46 


0.46 


Jun 


0.52 


0.55 


Jul 


0.51 


0.59 


Aug 


0.62 


0.66 


Sep 


0.80 


0.77 


Oct 


1.00 


0.99 



The upper bounds of the confidence intervals are 
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STX 


M3 


AT 

Nov 


1.42 


1.47 


Dec 


1.19 


1.77 


Jan 


0.79 


3.29 


Feb 


1.34 


2.32 


Mar 


0.64 


1.58 


Apr 


0.50 


0.55 


TV If 

May 


0.53 


0.53 


Jun 


0.60 


0.63 


T 1 

Jul 


0.59 


0.68 


Aug 


0.72 


0.76 


Sep 


0.93 


0.90 


Oct 


1.16 


1.15 



The statistics for the residual time series are 





mean 


stddev 


skewness 


kurtosis 


JB stats 


KS stats 


STX 


0.0000 


0.9987 


0.2997 


8.8683 


14306.7398 


0.1070 


M3 


0.0000 


0.9987 


0.8369 


12.9883 


31055.4814 


0.1213 



We see that the model fit is quite poor, which was expected with this long history of data. 

In the following table, we provide the correlation estimate and the corresponding lower and upper bounds 
of the confidence interval estimate: 





corr 


lower 


uppper 


Nov 


0.9244 


0.9077 


0.9383 


Dec 


0.6724 


0.6097 


0.7267 


Jan 


0.2439 


0.1458 


0.3372 


Feb 


0.7414 


0.6893 


0.7859 


Mar 


0.5302 


0.4529 


0.5995 


Apr 


0.8464 


0.8143 


0.8733 


May 


0.9585 


0.9493 


0.966 


Jun 


0.9456 


0.9335 


0.9555 


Jul 


0.9242 


0.9078 


0.9378 


Aug 


0.9401 


0.9271 


0.9509 


Sep 


0.9396 


0.9261 


0.9506 


Oct 


0.951 


0.9403 


0.9599 



We see that all the estimates are robust by looking at the confidence intervals. 
6.2 Spot-Prompt Model Daily Natural Gas Data 

Here we use around 2 year of data from 2/13/2008 to 11/25/2009 at two US locations STX and WLA together 
with the index data at those locations with the Spot-Prompt model. The data step is daily. We assume that 
the granularity of the term structures of mean-reversion level, volatility is monthly and correlation is fiat. 

The results of the calibration are summarized below. The mean-reversion rate estimates and their standard 
errors are given by 





STX 


WLA 


mr 


156.63 


163.99 


mrc 


19.79 


20.67 



We note that the mean-reversion rate estimates are typically higher when estimated with the Spot-Prompt 
model than with the one-factor mean-reversion model. 

We assume the step following interpolation of the term structures. The local volatility term structures are 
given by 
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STX 


STXIndcx 


WLA 


WLAIndcx 


Nov 


1.88 


0.82 


1.75 


0.74 


Dec 


0.93 


0.63 


0.79 


0.60 


Jan 


0.90 


0.51 


0.73 


0.50 


Feb 


0.72 


0.54 


0.70 


0.48 


Mar 


0.62 


0.81 


0.58 


0.69 


Apr 


0.47 


0.43 


0.48 


0.43 


May 


0.54 


0.70 


0.49 


0.70 


Jun 


0.56 


0.49 


0.57 


0.49 


Jul 


0.57 


0.72 


0.57 


0.74 


Aug 


0.67 


0.53 


0.67 


0.56 


Sep 


1.50 


1.28 


1.38 


1.31 


Oct 


1.82 


0.78 


1.77 


0.80 



The lower bounds of the confidence intervals are 





STX 


STXIndcx 


WLA 


WLAIndcx 


Nov 


1.53 


0.66 


1.42 


0.60 


Dec 


0.71 


0.48 


0.60 


0.46 


Jan 


0.68 


0.39 


0.56 


0.38 


Feb 


0.57 


0.43 


0.56 


0.38 


Mar 


0.51 


0.66 


0.48 


0.57 


Apr 


0.39 


0.36 


0.39 


0.36 


May 


0.44 


0.57 


0.40 


0.58 


Jun 


0.46 


0.40 


0.47 


0.41 


Jul 


0.47 


0.60 


0.47 


0.61 


Aug 


0.55 


0.44 


0.55 


0.46 


Sep 


1.23 


1.05 


1.13 


1.08 


Oct 


1.51 


0.64 


1.46 


0.66 



The upper bounds of the confidence intervals are 





STX 


STXIndcx 


WLA 


WLAIndcx 


Nov 


2.46 


1.07 


2.28 


0.97 


Dec 


1.37 


0.92 


1.16 


0.88 


Jan 


1.32 


0.74 


1.07 


0.74 


Feb 


0.95 


0.72 


0.93 


0.64 


Mar 


0.78 


1.02 


0.74 


0.88 


Apr 


0.60 


0.54 


0.60 


0.55 


May 


0.69 


0.89 


0.63 


0.90 


Jun 


0.71 


0.62 


0.72 


0.63 


Jul 


0.72 


0.91 


0.73 


0.93 


Aug 


0.85 


0.68 


0.86 


0.71 


Sep 


1.91 


1.63 


1.76 


1.67 


Oct 


2.30 


0.98 


2.23 


1.01 



The statistics for the residual time series are given by 





mean 


stddev 


skewncss 


kurtosis 


JB stats 


KS stats 


STX 


0.0191 


1.0018 


-0.4551 


1.9719 


88.8309 


0.0605 


STXIndcx 


-0.0126 


0.9876 


0.3218 


1.2814 


38.7256 


0.0503 


WLA 


0.0146 


1.0045 


-0.3954 


2.3440 


115.2543 


0.0565 


WLAIndcx 


-0.0151 


0.9876 


0.3836 


1.5635 


57.1260 


0.0583 



We see that the model fit is much better in this case than in the previous example. 
The correlation matrix is given by 
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STX 


STXIndcx 


WLA 


WLAIndcx 


1 a a 
1.1)1) 


A 1 O 

-O.lz 


a a a 
0.94 


A 1 O 

-O.lz 


-0.12 


1.00 


-0.14 


0.96 


0.94 


-0.14 


1.00 


-0.15 


-0.12 


0.96 


-0.15 


1.00 



and the corresponding lower bound of the confidence interval estimate 



STX 


STXIndcx 


WLA 


WLAIndcx 


1.00 


-0.21 


0.93 


-0.21 


-0.21 


1.00 


-0.23 


0.95 


0.93 


-0.23 


1.00 


-0.23 


-0.21 


0.95 


-0.23 


1.00 



and the upper bound of the confidence interval estimate 



STX 


STXIndcx 


WLA 


WLAIndcx 


1.00 


-0.03 


0.95 


-0.03 


-0.03 


1.00 


-0.05 


0.96 


0.95 


-0.05 


1.00 


-0.05 


-0.03 


0.96 


-0.05 


1.00 



6.3 Log-Normal Mean-Reversion Model Hourly Load and Power Data 

Here we use around 5.5 years of hourly ERCOT power locational marginal prices (LMP) and load data from 
7/1/2004 12am to 2/28/2010 11pm with the log-normal mean-reversion model. The data step is hourly. We 
assume that the granularity of the term structures of mean-reversion level, volatility, correlation is monthly. 
The valuation date is the next day after the the last date in the time series, i.e. 3/1/2010. 

The results of the calibration are summarized below. The mean-reversion rates and standard error estimates 
are given by 





Load 


LMP 


mr 


476.73 


668.39 


mre 


13.18 


15.86 



We assume the step following interpolation of the term structures. The local volatility term structures are 
given by 





Load 


LMP 


Mar 


8.37 


14.22 


Apr 


9.30 


14.15 


May 


9.32 


15.51 


Jun 


8.36 


15.17 


Jul 


7.91 


13.47 


Aug 


7.96 


13.73 


Sep 


8.94 


14.50 


Oct 


9.28 


15.89 


Nov 


8.78 


15.39 


Dec 


7.73 


14.76 


Jan 


7.31 


15.22 


Feb 


7.26 


14.19 



The lower bounds of the confidence intervals are 
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The upper bounds of the confidence i 





Load 


LMP 


Mar 


8.18 


13.91 


Apr 


9.09 


13.83 


May 


9.11 


15.17 


Jun 


8.17 


14.82 


T 1 

Jul 


7.75 


13.20 


Aug 


7.80 


13.45 


Sep 


8.76 


14.20 


Oct 


9.09 


15.57 


Nov 


8.60 


15.07 


Dec 


7.58 


14.46 


Jan 


7.16 


14.91 


Feb 


7.11 


13.89 


jrvals are 




Load 


LMP 


Mar 


8.56 


14.55 


Apr 


9.52 


14.49 


May 


9.54 


15.87 


Jun 


8.55 


15.53 


Jul 


8.08 


13.76 


Aug 


8.13 


14.02 


Sep 


9.14 


14.81 


Oct 


9.47 


16.23 


Nov 


8.97 


15.72 


Dec 


7.90 


15.08 


Jan 


7.46 


15.55 


Feb 


7.42 


14.51 



We note that the mean-reversion rates and volatilities are much higher for hourly power data than for daily 
gas data, this is a typical feature of the results of model calibration for power data. 
Here are the statistics for the residual time series 





mean 


stddev 


skewncss 


kurtosis 


JB stats 


KS stats 


Load 


0.0000 


0.9999 


0.0022 


-0.0532 


5.8915 


0.0116 


LMP 


0.0000 


0.9999 


0.7021 


3.1254 


24288.7583 


0.0673 



We see that the model fit is relatively good for the Load data. 

In the following table, we provide the correlation estimate and the corresponding lower and upper bounds 
of the confidence intcrnval estimate: 





corr 


lower 


uppper 


Mar 


0.706 


0.6895 


0.7218 


Apr 


0.638 


0.6182 


0.657 


May 


0.6142 


0.5938 


0.6338 


Jun 


0.7176 


0.7014 


0.7331 


Jul 


0.7294 


0.7153 


0.7428 


Aug 


0.7035 


0.6883 


0.718 


Sep 


0.649 


0.6313 


0.6659 


Oct 


0.653 


0.6359 


0.6695 


Nov 


0.6968 


0.6811 


0.7118 


Dec 


0.7233 


0.709 


0.737 


Jan 


0.6971 


0.6817 


0.7119 


Feb 


0.7006 


0.6846 


0.716 
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7 Conclusion 



We developed a simple robust approach for the joint historical calibration of several energy underlyings. This 
approach takes into account seasonality effects and uneven frequency of data. It allows to choose different 
granularity of model parameter term structures that would provide more robust estimates based on the com- 
puted confidence intervals of model parameters. We also provided a simple way to check the goodness of 
model fit. It shows whether the chosen models for each underlying are a good choice for the description of 
the underlying dynamics. A basic procedure of how to remove data outliers was also briefly mentioned. More 
elaborate methods for dealing with data outliers is an interesting question and can be a topic of a future study. 

We presented several examples of calibration results for several data sets. These examples provide a good 
illustration of typical parameter values and term structures for natural gas and power underlyings. This 
provides a good benchmark and guidance for calibration with other energy data sets. 
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